Ground state cooling of atoms in optical lattices 
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We propose two schemes for cooling bosonic and fermionic atoms that are trapped in a deep 
optical lattice. The first scheme is a quantum algorithm based on particle number filtering and state 
dependent lattice shifts. The second protocol alternates filtering with a redistribution of particles 
by means of quantum tunnelling. We provide a complete theoretical analysis of both schemes and 
characterize the cooling efficiency in terms of the entropy. Our schemes do not require addressing 
of single lattice sites and use a novel method, which is based on coherent laser control, to perform 
very fast filtering. 
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I. INTRODUCTION 

Ultracold atoms stored in optical lattices can be con- 
trolled and manipulated with a very high degree of pre- 
cision and flexibility. This places them among the most 
promising candidates for implementing quantum compu- 
tations Q, 0, S S an d quantum simulations of certain 
classes of quantum many-body systems [tl l?l IMI^. ITollTll 
However, both quantum simulation and quan- 
tum computation with this system face a crucial prob- 
lem: the temperature in current experiments is too high. 
In this paper we propose and analyze two methods to de- 
crease the temperature and thus to reach the conditions 
required to observe the interesting regimes in quantum 
simulations and quantum computation. 

So far, several experimental groups have been able to 
load bosonic or fermionic atoms in optical lattices and 

reach the strong interaction regime m e e m h 

IT9I I20I l2l| . In those experiments, the typical tempera- 
tures are still relatively high. For instance, the analysis 
of experiments in the Tonks gas regime indicates a tem- 
perature of the order of the width of the lowest Bloch 
band [16( , and for a Mott Insulator a temperature of the 
order of the on-site interaction energy has been reported 
[T7L 122^ . For fermions one observes temperatures of the 
order of the Fermi energy [2^, 0, ■ Those tempera- 
tures put strong restrictions on the physical phenomena 
that can be observed with those systems and also on the 
quantum information tasks that can be carried out with 
them. They stem from the fact that atoms are loaded adi- 
abatically starting from a Bose-Einstein condensate (in 
the case of bosons). On the one hand, the original con- 
densate has a relatively high entropy [26| that is inherited 
by the atoms in the lattice in the adiabatic process. On 
the other hand, the process may not be completely adi- 
abatic, which gives rise to heating. Thus, it seems that 
the only way of overcoming these problems is to cool the 
atoms once they have been loaded in the optical lattice. 

One may think of several ways of cooling atoms in 
optical lattices. For example, one may use sympathetic 
cooling with a different Bose-Einstein condensate fl3l 
|27^ . Here we propose two alternative schemes which do 



not require the addition of a condensate. They aim at 
cooling atoms to the ground state of the Mott-insulating 
(MI) regime and allow us to predict temperatures which 
are low enough for practical interests. Our protocols are 
based on translation invariant operations (i.e. do not 
require single-site addressing) and include the presence 
of an additional harmonic trapping potential, as it is the 
case in present experiments. Although we will be mostly 
analyzing their effects on bosonic atoms, they can also 
be used for fermions. 

The first scheme is based on the repeated application 
of occupation number filtering |2S[ . Via tunnelling, par- 
ticles from the borders of the trap are transferred to the 
center, where they are discarded by subsequent filter op- 
erations. The second scheme combines filtering with al- 
gorithms inspired by quantum computation J5| and hence 
will be termed algorithmic cooling of atoms [23 ] . The cen- 
tral idea is to split the atomic cloud into two components 
and to use particles at the border of one component as 
"pointers" that remove "hot" particles at the borders of 
the other component. We provide a detailed theoretical 
description of our cooling schemes and compare our the- 
ory with exact numerical calculations. In particular, we 
quantify the cooling efficiency analytically in terms of the 
initial and final entropy. We find that filtering becomes 
more efficient at low temperatures. This feature makes it 
possible to reach states very close to the ground state af- 
ter only a few subsequent filtering operations. Our theory 
further predicts that the algorithmic protocol is more ef- 
ficient at higher temperatures and that the final entropy 
per particle becomes zero in the thermodynamic limit. 
In addition, experimental requirements and time scales 
are discussed. 

Since filtering is an important ingredient of all our pro- 
tocols, we have devised an fast physical implementation 
which is based on optimal coherent laser control. Al- 
ready comparatively simple optimization schemes work 
on a time scale that is significantly shorter than the one 
in H3 and |3C| . 

The paper is organized as follows. We start in Sect. 
Illl with reviewing the physical system in terms of the 
Bose-Hubbard model and discussing realistic initial state 
variables such as entropy and particle number. In Sect. 
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IIIII we give a detailed theoretical analysis of filtering un- 
der realistic experimental conditions. Building on this, 
we study in Sect. HVI the repeated application of filtering. 
An algorithmic ground state cooling protocol is proposed 
and analyzed in Sect. [V] Next, Sect. IVII is dedicated to 
the discussion of our protocols, including a comparison of 
analytical results with numerical calculations. A further 
central result of our work is presented in Sect. IVIII There 
we propose an ultra-fast implementation of filtering op- 
erations based on coherent laser control. We conclude 
with some remarks concerning possible variants and ex- 
tensions of our protocols. In the Appendix we develop a 
fcrmionization procedure of the Bose-Hubbard Hamilto- 
nian, which accounts for up to two particles per site. 



II. INITIAL STATE AND BASIC CONCEPTS 

We consider a gas of ultra-cold bosonic atoms which 
have been loaded into a three dimensional (3D) optical 
lattice. The lattice depth is proportional to the dynamic 
atomic polarizability times the laser intensity. We further 
account for an additional harmonic trapping potential 
which either arises naturally from the Gaussian density 
profile of the laser beams or can be controlled se parately 
via an external magnetic or optical confinement |lq . |32| . 

In the following we will restrict ourselves to one- 
dimensional (ID) lattices, i.e. we assume that tunnelling 
is switched off for all times along the transversal lattice 
directions. This system is most conveniently described in 
terms of a single band Bose-Hubbard model [3j|. For a 
lattice of length L the Hamiltonian in second quantized 
form reads 



L/2 
k=-L/2 



t U 2 

-J(a\a k+l + h.c.) + —n k (n k - 1) + bk n k 



The parameter J denotes the hopping matrix element 
between two adjacent sites, U is the on-site interaction 
energy between two atoms and the energy b accounts for 
the strength of the harmonic confinement. Operators a\, 
and a/c create and annihilate, respectively, a particle on 
site k, and n k = a\.a k is the occupation number oper- 
ator. When raising the laser intensity the hopping rate 
decreases exponentially, whereas the interaction param- 
eter U stays almost constant j3^|. Therefore we have 
adopted U as the natural energy unit of the system. 

In the following we will consider ID thermal states in 
the grand canonical ensemble, which are characterized 
by two additional parameters, temperature kT = 1/(3 
and chemical potential p. In particular, we are inter- 
ested in the no-tunnelling limit |34j . J — > 0, in which the 
Hamiltonian Q becomes diagonal in the Fock basis of in- 
dependent lattice sites: {|«-(l-i)/2 • ■ ■ n o ■ ■ ■ n (L-i)/2)}- 
The density matrix then factorizes into a tensor product 



of thermal states for each lattice site: 

p=I e -<^-^)= ® p k , (2) 

fe=-L/2 

which simplifies calculations considerably. For instance, 
the von-Neumann entropy can be written as, 

(3) 



S{ P )=tx{p\og 2 p)=Y J S{p k ). 



This quantity will be central in this article, because it 
allows to assess the cooling performance of our proto- 
cols. To be more precise, we define two figures of merit. 
The ratio of the entropies per particle after and before 
invoking the protocol, Sf/st, quantifies the amount of 
cooling. The ratio of the final and initial number of par- 
ticles, Nf/Ni, quantifies the particle loss induced by the 
protocol. 

Note, however, that the entropy S(pf) is only a good 
figure of merit if the state pj after the cooling protocol 
is close to thermal equilibrium. If this is not fulfilled, 
we compute an effective thermal state, pf —* p e s, by ac- 
counting for particle number and energy conservation in 
closed, isolated systems. This is performed numerically 
by tuning the chemical potential and temperature of a 
thermal state p e g until the expectation values for par- 
ticle number and energy coincide with the ones of the 
original state pf. This procedure can be implemented 
rather easily in the no-tunnelling regime, in which the 
density matrix factorizes lf2"|. 

Our figures of merit can then be calculated from p c g. 
For instance, the final entropy is given by S(p e g). It 
constitutes the maximum entropy of a state which yields 
the same expectation values for energy E and particle 
number N as the final state pj. In this context it is im- 
portant to point out that other variables, like energy or 
temperature, are not very well suited as figures of merit, 
because they depend crucially on external system param- 
eters such as the harmonic trap strength. 

We now study the structure of the initial state in more 
detail. To this end we first give typical parameter values. 
The analysis of recent experiments in the MI regime im- 
plies a substantial temperature of the order of the on-site 
interaction energy 0, 0, \^ . This result is consistent 
with our own numerical calculations |40| and translates 
into an entropy per particle s := S/N ~ 1. The particle 
number in a ID tube of a 3D lattice as in 0] typically 
ranges between N = 10 and N = 130 particles. A repre- 
sentative density distribution corresponding to such ini- 
tial conditions (with N = 65) is plotted in Fig. In this 
example the inverse temperature is given by (3U = 4.5. 
Since our cooling protocols lead to even lower tempera- 
tures, we will from now on focus on the low temperature 
regime, (3U 3> 1. Moreover, we will only consider states 
with at most two particles per site, which puts the con- 
straint p < 2U — 1/(3 on the chemical potential. Such 
a situation can either be achieved by choosing the har- 
monic trap shallow enough or by applying an appropriate 
filtering operation pgj. 
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Under the assumptions e@ u ^> 1 and fi — U/2 > 
b + 1/(2/3) we will now show that the density distribu- 
tion of the initial state 101 can be separated into regions 
that are completely characterized by fermionic distribu- 
tion functions of the form: 



1 



I _|_ e l3(bk2-p,) 



(4) 



To be more precise, for sites at the borders of the 
density distribution, bk 2 > fi - U/2 + 1/(2/3), the 
mean occupation number is given by (rife) ~ ni(k) with 
ni(fc) := fk(b,/3,fi). In the center of the trap, bk 2 <C 
fi - U/2 - 1/(2(3), one has: (n k ) rs 1 + n n (k) with 
nu(k) := fk(b, 0, fin) and effective chemical potential 
flu := fi-U (see e.g. Fig. [Up- 
starting from state J2J, with parameters 0, u and b, 
the grand canonical partition function for site k is given 
by: 



e fe 



1 + a Xk + b x 2 , 



(5) 



with a;,, = e-^ fc , a = and 6 = e^ -17 ). In this 
notation the probabilities p£ of finding n particles at site 
k can be written as: pi = 1/Ofc, p\ — a x/Qk and 
p\ = b x 2 /Qk- For analyzing these functions we split 
the lattice into a central region and two border regions. 
For lattice sites at the borders one finds bx 2 -C 1, ax, 
meaning that the probability for doubly occupied sites 
becomes negligible: p\ <C p1-,p\- The average occupa- 
tion is thus given by (rife) « p\, with 



Pk 



1 



1 + ax 1 + e /3(bfc 2 -M) ' 



(6) 



In the crossover region, 6fc 2 i=s jti — C//2, one obtains a 
MI phase {pi, pi < ~ 1). In the center of the trap 
one finds a negligible probability for empty sites: pi <C 
p\,p\, since ax,bx 2 ^> 1. The average population at site 
k becomes (n k ) =pJ+2^Ril +Pk> where 



vl 



bx 2 



1 



ax - 



bx 2 



l + e i3{bk*-{ti-u)) ■ 



(7) 



This is identical to the fermionic distribution @ with 
effective chemical potential p — U . Hence the density 
distribution in this lattice region can be interpreted as 
a thermal distribution of hard-core bosons (phase II in 
Fig. |Ui) sitting on top of a MI phase with unit filling. 
Note that this central MI phase is well reproduced by 
the function ni(k) , which originally has been derived 
for the border region. As a consequence, the density 
distribution for the whole lattice can be put in the simple 
form: (nk) ~ ni(k) + nu(k), which corresponds to two 
fermionic phases I and II [Fig. QJi], sitting on top of each 
other. In other words, the initial state of our system 
can effectively be described in terms of non-interacting 
fcrmions, which can occupy two different energy bands 
I and II, with dispersion relations ei = bk 2 and Su = 
bk 2 + U, respectively [see also Appendix A and |43|]- 



The initial density profile can be further characterized 



by two distinctive points. At sites k — ±kn 



Ifi/b, 



which correspond to the Fermi levels of phase I, one ob- 
tains (rife ) = 1/2. Hence, k^ determines the radius of 
the atomic cloud. Note also that in the case fi w U 
singly occupied sites around the Fermi levels become de- 
generate with doubly occupied sites at the center of the 
trap. At the central site (k — 0) one one finds an average 
occupation: 



(n ) 



1 



1 



(8) 



For instance, the value (no) 
potential to be u = U. 



= 3/2 fixes the chemical 



III. ANALYSIS OF FILTERING 

By filtering we denote the state selective removal of 
atoms from the system, depending on the single site oc- 
cupation number [2^ . For instance, this can be achieved 
with a unitary operation 



TT M,m—M i n \ 



\M,m- M), 



(9) 



that transfers m — M particles from a Fock state with 
m atoms in internal states \a) to an initially unoccu- 
pied level 1 6). Particles in this level are removed subse- 
quently and the process is repeated for all m > M. The 
maximum single site occupation number then becomes 
M . Alternatively, filtering can be described in terms of 
a completely positive map acting solely on the density 
operator of atoms \a): 



F M ■ y^Pn,m|n)(r 



\M)(M\.(10) 

n,m<M n>M 

In particular, we are interested in the filtering operation 
F\ which yields either empty or singly occupied sites. 
This operation can be carried out with the scheme intro- 
duced in [2^] , which is based on the blockade mechanism 
due to atom-atom interactions. It enables a state se- 
lective adiabatic transfer of particles from one internal 
state of the atom to another. Recently, an alternative 
scheme which relies on resonant control of interaction 
driven spin oscillations has been put forward jsO]. How- 
ever, the predicted operation times of both approaches 
are comparatively long. Since fast filtering is crucial for 
the experimental realization of our protocols, we will pro- 
pose in Sect. IVIII an ultra-fast, coherent implementation 
of F\, relying on optimal laser control. 

In Fig. n we study the action of F\ on a realistic ID 
thermal state in the no-tunnelling regime, as defined in 
the previous section. One observes that a nearly perfect 
MI phase with filling factor v = 1 is created in the center 
of the trap. Defects in this phase are due to the presence 
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FIG. 1: Spatial dependence of the average particle number (n) 
and entropy S before (left) and after (right) the application 
of the filter operation Fx. The final particle distribution can 
be well described by Eq. JIJ (dashed line in (b)). Numerical 
parameters for initial state: Ni = 65, Sj=l and U/b = 700 
(f3U = 4.5, fi/U = 1). Figures of merit: Sf/si = 0.56 and 
N f /Ni = 0.80. 



of holes and preferably locate at the borders of the trap. 
This behavior is reminiscent of fcrmions, for which ex- 
citations can only be created within an energy range of 
order kT around the Fermi level. This numerical observa- 
tion can easily be understood with our previous analysis 
of the initial state. Filtering removes phase II, which is 
due to doubly occupied sites, and leaves the fermionic 
phase I unaffected [Fig. EKI- 

Let us now study the cooling efficiency of operation 
F\. This means we have to compute the entropy per 
particle of the states before and after filtering. According 
to our preceding discussion this problem reduces to the 
computation of the entropy S and the particle number 
N, corresponding to the bands I and II. The entropy of 
a fermionic distribution of the form Q is given by: 



S F (b,P,z) 



Li 3/2 (-2))] , (11) 



with fugacity z = e"^ and XA n {x) denoting the polylog- 
arithm functions. The function o{z) is defined as the 
integral 



a[z) 



dx log 



1 +z e 



(12) 



For phase I one can find a simpler expression for the 
entropy by expansion around the Fermi level k — 
k^ + dk. Note that the range of validity, \dk\ -C 1/s/^b, 
of this approximation covers all lattice sites that give a 
significant contribution to the total entropy. This yields 



the following relations: 



(13) 



with CI := 7r 2 /(61n2). For phase II this approach is typ- 
ically not valid and one obtains the general expressions: 



Sn = S F (b,f3, zn), N u = 



(14) 



with zn = e^ 11 



vim z\\ = In the special case fj, = U 

limplify the above expressions to: 

Su w cttt __ , Nn ~ Vn , ; 



one can 



(15) 



with numerical coefficients on sa 2.935 and rju = (1 — 

v^)\ZtFC(1/2) ~ 1-063. 

With these findings we can now give a quantitative 
interpretation of Fig. ^ The initial entropy is composed 
of two components: Si = Si + Sn. Filtering removes 
the contribution Sn, which arises from the coexistence 
of singly and doubly occupied sites. The final entropy is 
thus given by Sf = Si. This residual entropy is localized 



around the Fermi levels 
width A [Pig. El: 

A = 



-fc M and kfj, within a region of 



(16) 



and one can write St — <J\ A. For the initial and fi- 
nal particle numbers one has the corresponding relations: 
Ni = Ni + Nn and Nf = N\. Hence, the final entropy 
per particle can be written as: 



Sf. 

Nf 



CTl 



(17) 



For the special choice [i — U (or equivalently (no) = 1.5) 
one finds the following expressions for our figures of 
merit' 



it 

■Si 

N 



cri 



erf 

1 



mi 

2yW 



(18) 
(19) 



This result shows that filtering becomes more efficient 
with decreasing temperature, since Sf/si cx 1/ ' v / '/3U — * 
and Nf/N -> 1 for (3U -> oo. 

It is important to note that the state after filter- 
ing is not an equilibrium state, because it is energet- 
ically favorable that particles tunnel from the borders 
to the center of the trap, thereby forming doubly oc- 
cupied sites. According to the discussion in the pre- 
vious section this implies that the final entropy, which 
enters the cooling efficiency, should be calculated from 
an effective state p e // after equilibration. However, this 
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would yield a rather pessimistic estimate for the cool- 
ing efficiency. Since the final density nj(k) already has 
the form of a thermal distribution function, one can eas- 
ily come up with a much simpler (and faster) way to 
reach a thermodynamically stable configuration. While 
tunnelling is still suppressed, one has to decrease the 
strength of the harmonic confinement to a new value 
with b' < bU /(2/i). The system is then in the equilibrium 
configuration /&(&',/?', «') (@J with rescaled parameters 
V = T b'/b and fj,' = fi b'/b < U/2. This observation 
shows that it is misleading to infer the cooling efficiency 
solely from the ratio T'/T, because it depends crucially 
on the choice of b'. Note also that this procedure indeed 
allows to achieve the predicted value lfT%| for the cooling 
efficiency. 



IV. GROUND STATE COOLING WITH 
SEQUENTIAL FILTERING 

We have seen that filtering is limited by the fact that 
it cannot correct defects that arise from holes in a perfect 
MI phase. In order to circumvent this problem, we will 
now consider a repeated application of filtering, which 
will clearly profit from the increasing cooling efficiency 
as temperature is decreased. Our approach requires to 
iterate the following sequence of operations: (i) we allow 
for some tunnelling while the trap is adjusted adiabati- 
cally in order to reach a central occupation of (n ) w 1.5; 
(ii) we suppress tunnelling and perform the filtering op- 
eration Fi ; (iii) the trap is slightly opened so that the fi- 
nal distribution resembles a thermal distribution of hard- 
core bosons. This way we transfer "hot" particles from 
the borders to the center of the trap, where they can be 
removed by subsequent filtering. 

We are interested in the convergence of the entropy 
and temperature as a function of the number of itera- 
tions. However, the adiabatic process is very difficult to 
treat both analytically and numerically. Therefore we 
distinguish in the following between three different sce- 
narios that are based on specific assumptions. 

(i) Thermal equilibrium: We assume that the en- 
tropy is conserved and that the system stays in ther- 
mal equilibrium throughout the adiabatic process. Since 
this condition will in general not be fulfilled for closed, 
isolated quantum systems, the following analysis can 
only provide a rather rough description of the real sit- 
uation. To be more precise, we start from a thermal 
state with initial parameters (3, b and u. After filter- 
ing and adiabatic evolution one has a thermal state in 
the no-tunnelling regime with new parameters /?', b' and 
p! . As we have shown in the previous sections, ther- 
mal states in the no-tunnelling regime can effectively be 
described in terms of two fermionic components. This 
allows us to determine the new parameters f3' and b' 
by identifying: = Si(/3', b') + S n (/?',&') and 

Ni(J3,b) = Ni(0',b') + N u {f3',b'). The desired central 
filling (no) = 1.5 fixes the chemical potentials to be 



/i = // = U. Using expressions (|13(l and l|15(l one finds: 

P'U = [A+^A 2 + 4/3C/) 2 /4, (20) 

(21) 




with A = (<Tiij3U/<Ji — r]u)/2. After a second filtering 
operation the entropy per particle is thus given by: 



S2 = 0"! 



1 

Wu' 



(22) 



Since our analysis only holds in the limit f3U ^> 1, one can 
simplify the above expressions to: f3'U « (au(3U / (2a\)) 2 
and b' ss b. This allows us to establish a simple recursion 
relation for the entropy per particle s n after the n-th 
filter operation: 



4cri 



(23) 



Since the limit (3U ^S> 1 implies s < 1, one finds that the 
entropy per particle converges extremely fast to zero. 

(ii) Quantum evolution: Let us now study a more 
realistic situation. To this end we resort to an effective 
description of the quantum dynamics in terms of two 
coupled Fermi bands [Appendix A ]: 

H = [bk 2 c{c k + (bk 2 + U)d\d k 

k 

- J {c\c k+ i + h.c.) 

- V2J (cJ.dfe+1 + d\c k+1 + h.c.) 

- 2 J (dj.dfe+1 + h.c.)]. (24) 

Here, operators Cfc,c£ refer to energy band I and dk,dt 
to band II, which is shifted from the lower band by the 
amount of the interaction energy U (see Sect. [H] and 
Fig. |2| . This treatment is self-consistent as long as the 
probability of finding a particle-hole pair is negligible, i.e. 
{ckc\d\dk) ~ 0. In Sect. [H]we have shown analytically 
that this is indeed fulfilled for thermal states in the no- 
tunnelling regime and for low temperatures (j3U 1). 
We have checked that it also holds for thermal states 
at finite hopping rate J, provided that J is not too big 
(J < 0.5 U). Since the Hamiltonian (|24|l is quadratic, 
we can study the complete protocol in terms of a single- 
particle picture. This can easily be done, if one further 
assumes that no level crossings occur in the course of the 
adiabatic evolution. Then, the state at any time t can 
be computed by populating the single-particle energies 
of H(t) according to the initial probabilities (after fil- 
tering) in energetically increasing order. This method is 
illustrated in more detail in Fig. [21 

After the initial filtering step only states in the lowest 
energy band are occupied. The occupation probability 
is given by the fermionic distribution fk(b,f3,fi) Q. In- 
creasing the trap strength to an appropriate value b' > b 
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FIG. 2: Effective description of thermai states in the no- 
tunneliing iimit in terms of non-interacting fermions occupy- 
ing two energy bands. The dispersion relations are £i = bk 2 
and £n = bk 2 + U, where k denotes the lattice site and U is 
the interaction energy. Increasing the harmonic trap strength 
from b to b' increases the chemical potential to p! so that 
the population of the upper band becomes energetically fa- 
vorable. In the bosonic picture this process corresponds to 
the formation of doubly occupied sites. 

in the course of the adiabatic process makes it energeti- 
cally favorable to occupy also the second band. We find 
the state after returning to the no-tunnelling regime, p', 
by populating the energy levels, corresponding to the new 
trap strength b', in energetically increasing order accord- 
ing to the initial probabilities /&(&, 0, p)- At this point 
we distinguish between two further scenarios: (ii.a) The 
state p' is mapped to a thermal state in the usual way 
by accounting for energy and particle number conserva- 
tion [Sect. HI) . This way we can quantify the amount 
of "heating" resulting from the fact that the system is 
not in thermal equilibrium at the end of the adiabatic 
process due to the different structure of the energy spec- 
trum, (ii.b) The next filtering operation acts directly on 
the time evolved state p 1 . 




iterations iterations 



tocol. In scenario (ii.b) the final entropy saturates at a 
finite value and the system is not in perfect thermal equi- 
librium. However, the final state still resembles a thermal 
state of hard-core bosons in a harmonic trap. 

These results imply that sequential filtering can clearly 
profit from equilibration. The reason is that equilibration 
reduces the defect probability in the center of the lower 
band and transfers entropy to the upper band, where it 
can be removed subsequently. This process in combina- 
tion with the increasingly high cooling efficiency of filter- 
ing at low temperatures can easily compensate the heat- 
ing induced by the adiabatic quantum evolution. From 
our data we can deduce that this heating corresponds 
to an entropy increase of around 20 % 37]. Without 
equilibration sequential filtering becomes very inefficient 
after the fourth iteration, which is also reflected in the 
excessive particle loss [Fig. |3J. The minimal attainable 
entropy is determined by the initial defect (hole) prob- 
ability in the center of the trap. Starting from a much 
colder state, which exhibits almost unit filling in the cen- 
ter of the lower band, therefore yields a final state very 
close to the ground state [Fig. 0. 

Remember that scenario (ii.b) is based on the assump- 
tion that no level crossings occur during the evolution. 
From our numerical analysis of the energy spectrum of 
(|24[1 we know, however, that level crossings can indeed 
appear (see also 0, H3). The reason is the vanishing 
small spatial overlap between single particle states at the 
border of the lower band and the center of the upper 
band. This has the following consequences for our previ- 
ous analysis: For rather small particle numbers (N < 15) 
level crossings are rare and inter-band coupling occurs al- 
ready for hopping rates, which are well within the range 
of validity of our single-particle description. We there- 
fore expect that our predictions, as depicted in Fig. |3 
are reasonable. For larger systems one has to tune the 
tunnelling rate deep into the supcrfluid regime J > 0.5 U 
in order to couple the two bands and to form doubly oc- 
cupied sites. However, this regime is no longer accessible 
within our fermionic model l|24|) . It remains to be inves- 
tigated to what extent this will alter our predictions for 
the cooling efficiency of sequential filtering. 



FIG. 3: (Color online) Entropy per particle S/N (left) and 
normalized number of particles N/No (right) as a function 
of the number of filtering iterations for initially iVo = 100 
particles. As discussed in the text, we distinguish between 
the scenarios (i) (black line), (ii.a) (red line) and (ii.b) (blue 
and cyan). 

Our results for all three cases are summarized in Fig. [3] 
We have computed numerically exact the entropy per 
particle as a function of the number of filtering cycles. 
Starting with an initial entropy sq — 1, scenarios (i) and 
(ii.a) predict that an entropy value close to zero can be 
obtained after only four iterations of the protocol [45| . 
According to our underlying assumptions the system is 
in thermal equilibrium after each iteration of the pro- 



V. ALGORITHMIC GROUND STATE COOLING 

A. The protocol 

We now propose a second cooling scheme, which we call 
algorithmic cooling, because it is inspired by quantum 
computation. As before the goal is to remove high energy 
excitations at the borders of the atomic cloud, which have 
been left after filtering. In contrast to sequential filtering 
we now restrict ourselves to a sequence of quantum op- 
erations that operate solely in the no-tunnelling regime. 
The central idea is to make use of spin-dependent lattices. 
A part of the atomic cloud can then act as a "pointer" 
in order to address lattice sites which contain "hot" par- 
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tides. In this sense the scheme is similar to evaporative 
cooling, with the difference that an atomic cloud takes 
the role of the rf-knife. Another remarkable feature of the 
protocol is that the pointer is very inaccurate in the be- 
ginning (due to some inherent translational uncertainty 
in the system), but becomes sharper and sharper in the 
course of the protocol. 
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FIG. 4: (Color online) Illustration of the protocol. The state 
is initialized with the filter operation Fi. (a) Particles from 
doubly occupied sites are transferred to state \b) (red) using 
operation ?7 2 'q. (b) The lattice |6) has been shifted 2k e sites 
to the right, so that the two distributions barely overlap, (c) 
Density distribution after k a lattice shifts. After each shift 
doubly occupied sites have been emptied. Afterwards lattice 
1 6) is shifted 4k e — k s sites to the left and an analogous filter 
sequence is applied, (d) The final distribution of atoms in 
state | a) is sharper compared to the initial distribution (dot- 
ted). Numerical parameters: Ni — 65, Si = 1, U/b = 300, 
fc e =21, k s = 20, N f = 30.2, s f = 0.31 (after equilibration). 

The steps of this algorithmic protocol are: (i) We begin 
with a thermal equilibrium cloud with two or less atoms 
per site, all in internal state |a), and without hopping. 
This can be ensured with a filtering operation F%. (ii) 
We next split the particle distribution into two, with an 
operation f7 20 [Fig-EKI- ( m ) The two clouds are shifted 
away from each other until they barely overlap. Then we 
begin moving the clouds one against each other, emptying 
in this process all doubly occupied sites. This sequence 
sharpens the density distribution of both clouds. It is 
iterated for a small number of steps, of order A i|16|) 
[Fig.^J;]. (iv) The atoms of type \b) are moved again to 
the other side of the lattice and a process similar to (iii) 
is repeated [Fig. 0JI] . (v) Remaining atoms in state |6) 
can now be removed. 

The final particle distribution cannot be made arbi- 
trarily sharp [Fig.^Ji], due to the particle number uncer- 
tainty in the tails of the distribution. In the following we 
will consider this argument more rigorously and develop 
a theoretical description of the protocol. 



B. Theoretical description 

For the sake of simplicity we consider a slightly mod- 
ified version of the protocol. The particle distributions 
| a) and \b) are now two identical but independent dis- 
tributions of hard-core bosons of the form (QJ [33 ■ The 
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FIG. 5: (Color online) Schematic description of the initial 
state for lattice sites k > 0: Two identical density distribu- 
tions for hard-core bosons, belonging to different species A 
(black) and B (blue), are shifted 2k e lattice sites apart. The 
region of non-integer filling has width 8 e . In the course of the 
protocol the lattice of species B is shifted k 3 = 35 e sites to 
the left. 

lattice \b) is shifted 2k E sites to the right. For given e 
the value k E = J Jijb^fl — In e//3p, is chosen such that 
for atoms |a) it holds (n£ } = e. This initial situation 
is depicted schematically in Fig. [5] The cutoff e defines 
also the width of the region with non-integer filling: 

5 E = yJpTJb (-/I - lne//3/i - a/I + lne//?/i) . (25) 

We analyze first a protocol that involves k s = 3<5 e lat- 
tice shifts and after each shift doubly occupied sites are 
emptied. Our goal is to compute the final shape of the 
density profile of atoms in state |a) (red line in Fig. [SJ. 
It is sufficient to consider only the reduced density ma- 
trices p a and p , which cover the range k £ [k £ — 2S E ; k £ ] 
and k G [k E ;k E + 2S E ], respectively. These density matri- 
ces can be written in terms of convex sums over particle 
number subspaces: 

26 e 

Pa = ]T Pa(Na)pa(N a ), (26) 
JV o =0 
25 e 

Pb = Pb(N b )p b (N b ). (27) 

^6 = 

The further discussion is based on the following cen- 
tral observation. If a state p a (N a ) interacts with a state 
Pa{N ) then our protocol produces a perfect MI state 
p'a(N a ) composed of N a = (N a ~N b )/2 particles. The fac- 
tor 1/2 arises from the fact that k s lattice shifts remove 
at most k s /2 particles from distribution \a). Note that 
this relation also allows for negative particle numbers, 
because N a merely counts the number of particles on the 
right hand side of the reference point k = k E — i/28 E . 
The final density matrix after tracing out particles in |6) 
can then be written as a convex sum over nearly perfect 
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(up to the cutoff error e) MI states 



S e /2 



P'a= E PaiKKiK), (28) 
N' a =-S e /2 



with probabilities 



Pa(K) 



2S, 

; E 

N b =6 e 



p a {2N' a + N b )p b {N b ). (29) 



The factor two is due to the fact that states with 
N a - N b = 2AI and N a - N b = 2AI + 1 are collapsed 
on the same MI state with N' a = M. Since Lyapounov's 
condition |4lj | holds in our system, we can make use of 
the generalized central limit theorem and approximate 
p a (N) = p b {N) by a Gaussian distribution with vari- 
ance a 2 = SN 2 = A/4 = 1/(2/375/7). Evaluation of 
Eq. (|29H then yields a Gaussian distribution with vari- 
ance a' 2 = a 2 1 '2. Since MI states do not contain holes, 



one can infer the final density distribution 



directly 



from p' a (N) by simple integration. This distribution can 
then be approximated by the (linearized) thermal distri- 
bution: 



1 



{n k ) 



1 + e 4(fc-fco)/ 



A' 



(30) 



The new effective tail width A' = V / A~7r/2 of the dis- 
tribution is roughly the square root of the original 
width A (|16|l . This effect leads to cooling, which we will 
now quantify in terms of the entropy. 

When applying similar reasoning also to the left side of 
distribution |a) one ends up with a mixture of MI states, 
which differ by their length and lateral position. This 
results in an extremely small entropy of the order Smi ~ 
log 2 A. However, this final state is typically far from 
thermal equilibrium. In order to account for a possible 
increase of entropy by equilibration, we have to compute 
the entropy of a thermal state, which has the same energy 
and particle number expectation values as the final state. 
In our case this is equivalent to computing the entropy 
directly from the density distribution l|3U|) : 



S' » at A' = 



aiy/n 



1 



V2 (/3 2 ^) 1/4 



N 



S. 



(31) 



where N = Nj and S — Si (|13|l correspond to the 
expectation values after the initial filtering operation. 
The final particle number is given by: N' ~ 2fco ~ 
N(l + \ne/(3fi). 

A significant improvement can be made by shifting the 
clouds only k s = 28 E sites. This prevents inefficient parti- 
cle loss, which has been included in our previous analysis 
in order make the treatment exact. With this variant the 
final particle number increases to N" = 2 (k e — S £ ), while 
in good approximation the final entropy is still given by 



S' . Hence, the final entropy per particle can be lowered 
to: 



S' 



N'' 



/7T 

2 I 



f 



(32) 



In s 

2/3m 



N 



This expression, which holds strictly only in the limit 
j3U ^> 1, shows that for fixed N the ratio s'/s becomes 
smaller at higher temperatures. Even more important, 
for fixed (3U, the entropy per particle s' decreases with 
as the initial number of particles in the system 
increases. 

Finally, let us remark that the final entropy can be 
further reduced, when the protocol is repeated with two 
independent states of the form (|28|) . In practice, this 
could be achieved with an ensemble of non-interacting 
atomic species in different internal levels. According to 
Eq. 129fl . each further iteration of the protocol decreases 
the total entropy by a factor 1/V2. 



VI. DISCUSSION OF RESULTS 

Let us now discuss and compare our previous results. 
In particular, we are interested in checking the range of 
validity of our analytical results by comparison with ex- 
act numerical calculations. To this end, we fix two pa- 
rameters, b/U and [i/U, and compute the entropy per 
particle as a function of the inverse temperature (3U 
[Fig. . 

Our results can be summarized as follows. Firstly, we 
find that our theoretical description is very accurate in 
the low temperature limit [3U 3> 1, and even holds in the 
(relatively) high temperature range f3U > 1. Secondly, 
our algorithmic protocol outperforms filtering consider- 
ably, especially in the high temperature range. Finally, 
based on the assumptions that underlie our calculations, 
subsequent filtering is typically superior to algorithmic 
cooling with respect to the minimal achievable entropy. 

We now discuss advantages, experimental require- 
ments and time scales of our cooling protocols. 



A. Sequential Altering 

(i) Advantages: If one combines filtering with equili- 
bration then the entropy per particle converges very fast 
to zero with the number of filter steps. The minimal 
value is only limited by finite size effects. Without equi- 
libration, the minimal entropy is limited by the finite 
probability of finding a hole in the central MI phase of 
the initial distribution. Furthermore, sequential filter- 
ing naturally allows for cooling in a 3D setup, because 
it preserves spherical symmetry. Note that filtering, and 
hence sequential filtering, can also be applied to fermionic 
atoms in an optical lattice |28|. 

(ii) Requirements and limitations: The repeated creation 
of doubly occupied sites in the center of the trap requires 
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FIG. 6: (Color online) Analytical (solid) and numerical (dots) 
results for the entropy per particle s = S/N versus the inverse 
temperature f3U for fixed trap strength b/U = 1/2500 and 
fixed fj, = U. We plot the initial value (black line) and the 
values after filtering with F\ (Eq. $17\ . brown line). This is 
compared to more elaborate cooling protocols: entropy after 
two iterations of sequential filtering including equilibration 
(Eq. 1221 , red line) , minimal entropy after sequential filtering 
without equilibration (cyan), and entropy after algorithmic 
cooling (Eq. I32L blue line). For the algorithmic protocol 
we have chosen e = 0.03 and k s = 2 S e . Note that this pro- 
tocol creates classical correlations between lattice sites. The 
numerics are therefore based on a representation of classical 
density matrices in terms of matrix-product states |4(i|l . 



precise control of the harmonic confinement over a wide 
range of values. In addition, non-adiabatic changes of lat- 
tice and external potentials might induce heating, which 
could reduce the cooling efficiency considerably, 
(iii) Time scales: The limiting factor here is not filtering 
but the adiabaticity criterion for changes of the hopping 
rate and the harmonic confinement. We have shown that 
after filtering the density distribution can be identified 
with a thermal distribution of spin-less fermions. Hence 
it is possible to find estimates for adiabatic evolution 
times based on single particle eigenf unctions, as calcu- 
lated in [3]]]. In the MI regime particle transfer from 
the borders to the center of the trap is very unlikely, be- 
cause the eigenfunctions of the upper and lower Fermi 
band do not overlap. Therefore, we propose as a first 
step to decrease the lattice potential at fixed harmonic 



confinement until eigenfunctions start overlapping. This 
process can still be treated within a fcrmionic (or Tonks 
gas) picture, since only the lower band is populated [Fig. 
12] • The average energy spacing around the Fermi level is 
SE = bN 2 U/N in the no-tunnelling regime, and 
stays roughly constant when passing over to the tun- 
nelling regime [3l|. As a consequence, the lattice po- 
tential should be varied on a time scale T ^> H/8E sa 10 
ms for N = 50 and h/U = 396 us |35|. For the sec- 
ond process, which involves the change of the harmonic 
potential at fixed hopping rate, it is more difficult to 
make analytic predictions for adiabatic time scales, since 
our single-particle description flty is no longer justified 
in general. One can, however, obtain a lower bound 
by considering the energy spectrum after returning to 
the no-tunnelling regime [Fig. [2J. The average energy 
spacing around the Fermi level is now dominated by 
the energy spacing at the bottom of the upper band: 
IW = vW/3/2 w U/(Ny/pU). This implies adiabatic 
evolution times which are a factor y / f3U/2 larger than 
for the first process. We have also verified the whole 
process numerically, using the matrix-product state rep- 
resentation of mixed states For initially N = 11 
particles we find adiabatic evolution times of the order 
T ~ 30 h/U for the first process, which is consistent with 
our analytical estimate. The second process is more time 
consuming with T > 120 H/U. 



B. Algorithmic cooling 

(i) Advantages: The algorithmic protocol operates 
solely in the no-tunnelling regime. Adiabatic changes 
of the lattice and/or the harmonic potential, which are 
time consuming and might induce heating, are therefore 
not required. Moreover it does not demand arbitrary 
control over the harmonic confinement. The correct ini- 
tial conditions can always be generated by the filter op- 
eration F\. The protocol is more efficient in the high 
temperature range and for large particle numbers. Ad- 
ditionally to ground state cooling, the algorithm can be 
used to generate an ensemble of nearly perfect quantum 
registers for quantum computation. This state, which 
can also be considered as an ensemble of possible ground 
states in the uniform system, might already be sufficient 
for quantum simulation of certain spin Hamiltonians. Fi- 
nally note that this protocol can naturally be applied also 
to fermionic systems. 

(ii) Requirements and limitations: The heart of the pro- 
tocol is the existence and control of spin-dependent lat- 
tices. Moreover, the algorithm is explicitly designed 
for cooling in one spatial dimension. Generalizations to 
higher dimensions are possible, but will typically not pre- 
serve the spherical symmetry of the initial density distri- 
bution. Moreover, one should keep in mind that the final 
states are typically far from thermal equilibrium. 

(iii) Time scales: Adiabatic lattice shifts can be per- 
formed very fast on a time scale determined by the on- 
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site trapping frequency v ~ 30 kHz. The limiting factor 
is the number of filter operations, which is of the or- 
der S £ (125(1 . Under realistic conditions this can amount 
to 50 operations. Filter operation times based on the 
adiabatic scheme l28| are of the order Tp ~ 200 h/U. 
With h/U — 396 us !,:j5j one finds a total operation time 
T ~ 630 ms, which is comparable with the typical par- 
ticle life time in present setups using spin-dependent lat- 
tices Q. We have studied this problem with an alterna- 
tive implementation of filtering [Sect. I VII) . which allows 
to reduce operation times by a factor of ~ 15, and hence 
makes algorithmic cooling feasible in current experimen- 
tal setups. 

VII. ULTRA-FAST FILTERING SCHEME 

We now propose an ultra- fast, coherent implementa- 
tion of filtering, which is based on optimal laser control. 
We restrict our discussion to the filtering operation Fi 
which is most relevant for our cooling protocols discussed 
above. We consider atoms in a particular internal level, 
| a), which are coupled to a second internal level, \b), via 
a Raman transition with Rabi frequency f2(t). In con- 
trast to the adiabatic scheme [2^ we consider constant 
detuning, but vary Q(t) in time. The Hamiltonian for a 
single lattice site reads 

H = ^-n a (h a - 1) + U ab h a h b + ^-n b (n b - 1) 

-(f2(t)a t 5 + fl*(*)S t o), (33) 

where U a , U b and U ab denote the on-site interaction ener- 
gies, according to the different internal states. Note that 
f2(t) can be complex, thus allowing for time-dependent 
phases. Our goal is to populate state \a) with exactly 
one particle per site which can be expressed by the 
unitary operation U : \N,0) -> \1,N - 1), V N e 
{1, 2, . . . , N max }. In order to do this, we control the 
time-dependence of fi(i) coherently and in an optimal 
way. To be more precise, we optimize a sequence of AI 
rectangular shaped pulses of equal length: 

M 

i=i 

After time T the system has evolved according to the 
unitary operator U(T). We want to minimize the de- 
viation of U(T) from the desired operation Uq, which 

we quantify by the infidelity e(T) = Y1n=T £n > wnere 
e N = 1 - |(1, N - 1\U(T)\N, 0)|. Since we allow for com- 
plex Rabi frequencies, e(T) is a function of 2M parame- 
ters {fi;, f2*} with I = 1, . . . , M. For given M and time T 
we optimize the cost function e(T) numerically using the 
Quasi-Newton method with a mixed quadratic and cubic 
line search procedure. This is repeated for different times 
T, while keeping the number of pulses M constant. We 
then increase M and repeat the whole procedure in order 
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FIG. 7: Operation error e vs. time tU a for different interaction 
energies: U b = 0.2U a , U ab = 0.2U a , S = 0.8 (solid); U b = U a , 
U ab = 0.6U a , 5 = 0.8 (dashed); U b = U a , U ab = 0.2U a , 
S = 1.6 (dotted); U b = U a = U ab , 5 = (dashed-dotted). 
We optimize a sequence of M = 10 pulses H34H . 

to check convergence of our results. In Fig.[7|we plot the 
minimal error e(T) for different interaction strengths. Al- 
ready for our simple control scheme we obtain very small 
errors, e ~ 10~ 4 , for a time T w 7/U a and interactions 
U a b = U b — Q.2U a . In comparison, the adiabatic scheme 
|28| would require T f» 150 /U a for the same set of pa- 
rameters. However, in deep contrast to (2Sj our method 
works for a very broad range of interaction energies and 
in addition allows for very high fidelities. 

It is important to remark that the operation time in- 
creases for small interaction anisotropics 6 = \U a + U b — 
2U ab \/U a and for 6 = our method fails. In the spe- 
cial case V ' a — U b = U ab this follows from the fact that 
in the Hamiltonian ilM-SI) the interaction part commutes 
with the coupling part. These problems can be solved, 
either by displacing the lattices that trap atoms \a) and 
1 6) and thereby reducing the effective interaction, U ab , or 
by performing more elaborate controls than the one from 
Eq. (J33J. 



VIII. CONCLUSION 

We have given a detailed analytical analysis of filter- 
ing in the no-tunnelling regime and in the presence of a 
harmonic trap. We have found that the residual entropy 
after filtering is localized at the borders of the trap, quite 
similar as in fermionic systems. Inspired by this result 
we have proposed two protocols that aim at removing 
particles from the borders and thus lead to cooling. One 
scheme transfers particles from the borders to the center 
of the trap, where they can be removed by subsequent fil- 
tering. In the other, algorithmic protocol particles from 
the system itself take the role of the rf-knife in evapora- 
tive cooling and remove directly particles at the borders. 
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We have quantified the cooling efficiency of these proto- 
cols analytically in terms of the initial and final entropy 
per particle. To this end, we have also developed an effec- 
tive description of the single-band Bose- Hubbard model 
in terms of two species of non-interacting fermions. 

A special virtue of our schemes is that they rely on 
general concepts which can easily be adapted to differ- 
ent experimental situations. For instance, our protocols 
can easily be extended to fermionic systems (for details 
see H3). Moreover, the algorithmic protocol can be im- 
proved considerably by the use of an ensemble of non- 
interacting atomic species in different internal states. In 
this context one should also keep in mind that a 3D lat- 
tice structure offers a large variety of possibilities, which 
have not been fully explored yet. 

Since we have identified the limitations of filtering, one 
can immediately think of alternative or supporting cool- 
ing schemes. For instance, one might use ring shaped 
lasers beams to remove particles at the borders of a 3D 
lattice. Equivalently, one can use a transition, which is 
resonant only for atoms with appropriate potential en- 
ergy [iiij. Although these methods do not allow to ad- 
dress individual lattice sites, they might still be useful as 
preliminary cooling steps for the protocols proposed in 
this article. 

We believe that the ideas introduced in this article 
greatly enhance the potential of optical lattice setups for 
future applications and might pave the way to the exper- 
imental realization of quantum simulation and adiabatic 
quantum computation in this system. We also hope that 
our analytical analysis of the virtues and limitations of 
current proposals, especially filtering, might trigger fur- 
ther research in the direction of ground state cooling in 
optical lattices. 
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APPENDIX A: MAPPING TO TWO-BAND 
FERMIONIC MODEL 

We start from the single-band Bose-Hubbard Hamil- 
tonian Qy and restrict the occupation numbers at each 



lattice site k to n k €E {0, 1,2}. In this truncated basis the 
Hamiltonian reads: 

H = J2 [bk 2 \l) k (l\ + 2(bk 2 + U)\2) k (2\ 

k 

- J (|0) fc |l) fe+1 (l| fe (0| fc+1 +h.c.) 

- v/2J (|0) fc |2) fc+ i<l|fc(l| fc+ i+h.c.) 

- V2J (|2) fe |Q) fe+1 (l| fe (l| fe+1 +h.c.) 

- 2J (|2) fc |l) fe+ i(l| fe (2| fe+ i+h.c.)]. (Al) 

One can now embed the three dimensional single site 
Hilbcrt space Hb = C 3 into the composite Hilbert space 
Hf = C 2 <g> C 2 of two species of hard-core bosons by 
applying the following mapping: 

|2) = ^(a t ) 2 |vac) c^vac), 
v2 

|1) = a^|vac) — > c^vac). 

(A2) 

Note that singly occupied bosonic states are mapped ex- 
clusively to the c-manifold, i.e. we omit the possibility 
of having one particle in the the d-manifold and no par- 
ticle in the c-manifold on the same site. After trans- 
forming hard-core bosons to fermions, c, d — > c, d, via a 
Jordan- Wigner transformation one obtains the following 
fermionic Hamiltonian: 

H = [bk 2 c\c k d k dl + {bk 2 + U)c\c k d\d k 

k 

- J (& k Ck+\dkd\dk+xd\ +1 +h.c.) 

- V2J (c\d k+ id k d) k c\ +l c k+ i+}i.c.) 

- V2J (d\c k+ ic k c\d k+ id\ +l +Yi.c.) 

- 2J (d' k dk+iCkc' k Ck+ic' k+1 + h.c.)] (A3) 

This Hamiltonian can also be written in the form H = 
P^HP, where H is the quadratic Hamiltonian l|24|) 
and P denotes the projection on the subspace, which 
is defined by c\d k = 0. This implies that bosonic 
atoms in an optical lattice can effectively be described 
in terms of the quadratic Hamiltonian Q24|). given that 
the probability of finding a particle-hole pair is negligi- 
ble, i.e. {c k c\d\d k ) w 0. 

Let us point out that similar fermionization procedures 
have been discussed in 0> El 
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